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Probability-Changing Cluster Algorithm: Study of Three-Dimensional 
Ising Model and Percolation Problem 
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We present a detailed description of the idea and procedure for the newly proposed Monte 
Carlo algorithm of tuning the critical point automatically, which is called the probability- 
changing cluster (PCC) algorithm [Y. Tomita and Y. Okabe, Phys. Rev. Lett. 86 (2001) 572]. 
Using the PCC algorithm, we investigate the three-dimensional Ising model and the bond 
percolation problem. We employ a refined finite-size scaling analysis to make estimates of 
critical point and exponents. With much less efforts, we obtain the results which are consistent 
with the previous calculations. We argue several directions for the application of the PCC 
algorithm. 
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1. Introduction 

The Ising model 1 ) is a basic model for studying phase 
, transitions and critical phenomena. The q-state Potts 

model, 2,3 ) which has q components for the order param- 
, eter, is a generalization of the Ising model. Then, the 
' Ising model corresponds to the q = 2 Potts model. The 

percolation model 4 ) also exhibits a geometric phase tran- 

■ sition. The Kasteleyn-Fortuin (KF) cluster representa- 
tion 5 ) of the q-state Potts model bridges the Potts model 
and the percolation problem; the bond percolation prob- 

, lem can be regarded as the q — 1 Potts model. Recently, 
' based on the cluster formalism, the multiple-percolating 
clusters of the Ising system with large aspect ratio have 

■ been studied. 6 ) 

For two-dimensional (2D) systems, exact or rigorous 
results have been obtained for the Potts models, 7, 8 ) and 
, they are used as the testing ground for numerical study. 

■ On the other hand, for three-dimensional (3D) systems, 
it is rare that exact results are available, and we rely 
more on numerical studies for revealing the nature of the 
problem. The Monte Carlo simulation 9 ) is a standard 
powerful tool to study critical phenomena numerically. 
To obtain accurate data, the development of efficient al- 

■ gorithms is highly demanded. Cluster algorithms 10 ' 11 ) 
, are examples of such efforts, and they have been suc- 

■ cessfully used to overcome slow dynamics in the Monte 
Carlo simulation. Swendsen and Wang (SW) 10 ) applied 

■ the KF 5 ) representation to identify clusters of spins. 

Extending the SW algorithm, we have recently pro- 
posed an effective cluster algorithm of tuning the crit- 
ical point automatically; this algorithm is called the 
probability-changing cluster (PCC) algorithm. 12 ) We 
have shown the effectiveness of the PCC algorithm for 
the case of 2D Potts models in the Letter. 12 ) The basic 
idea of our algorithm is that we change the probabil- 
ity of connecting parallel spins p in the KF representa- 
tion during the process of the Monte Carlo spin update. 
We decrease or increase p depending on the observation 
whether the KF clusters are percolating or not percolat- 
ing; essentially, we change the temperature. This simple 



negative feedback mechanism together with the finite- 
size scaling (FSS) 13 ) property of the existence probability 
(also called the crossing probability) E p , the probability 
that the system percolates, leads to the determination of 
the critical point. Since our ensemble is asymptotically 
canonical as Ap, the amount of the change of p, becomes 
0, the distribution functions of physical quantities obey 
the FSS; as a result, we can determine critical exponents 
using the FSS analysis. 

Previously, Machta et a/. 14 ) proposed another idea 
of cluster algorithm to tune the critical point automati- 
cally, which is called the invaded cluster (IC) algorithm. 
However, the ensemble of the IC algorithm is not neces- 
sarily clear, and it has a problem of "bottlenecks" , which 
causes the broad tail in the distribution of the fraction 
of the accepted satisfied bonds. 14 ) In contrast, it is guar- 
anteed that we approach the canonical ensemble in our 
PCC algorithm. 

In this paper we give the more detailed description of 
the PCC algorithm, and show the results for the 3D Ising 
model and the 3D bond percolation problem. We pay 
attention to the refined FSS analysis for determining the 
critical point and critical exponents, which is the same 
idea as was used in a high-resolution Monte Carlo study 
by Ferrenberg and Landau. 17 ) The rest of this paper is 
organized as follows: In §2 we explain the idea and the 
procedure of the PCC algorithm. In §3 the results for 
the 3D Ising model and percolation problem are shown, 
and the refined FSS analysis is discussed. Finally, we 
summarize this paper and give discussions in §4. 

2. Probability-Changing Cluster Algorithm 

2. 1 Idea of PCC algorithm for Ising model 

We start with explaining the idea of the PCC algo- 
rithm. Here, we deal with the ferromagnetic Ising model, 
whose Hamiltonian is given by 
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where J is the exchange coupling constant, and the sum- 
mation is taken over the nearest-neighbor pairs (i, j). In 
this case, the probability of connecting parallel spins in 
the KF representation is given by p = 1 — e _2J / fcBT . The 
procedure of Monte Carlo spin update is as follows: 

1. Start from some spin configuration and some value 
of p. 

2. Make KF clusters using the probability p, and check 
whether the system is percolating or not. Update 
spins following the same rule as the SW algorithm, 
that is, flip all the spins on any KF cluster to one 
of two states. 

3. If the system is percolating (not percolating) in the 
previous test, decrease (increase) p by Ap (> 0). 
Essentially, we change the temperature T. 

4. Go back to the process 2. 

Since we use the cluster representation and assign clus- 
ters, we are ready to check whether the system is per- 
colating or not in the process 2. The distribution of p 
for Monte Carlo samples approaches the Gaussian dis- 
tribution of which mean value is p c (L), after repeating 
the above processes. Here, p c {L) is the probability of 
connecting spins, such that the existence probability E p 
becomes 1/2, and depends on the system size L. The 
existence probability E p is the probability that the sys- 
tem percolates. In the limit of Ap — > 0, we approach the 
canonical ensemble, which will be discussed later. Then, 
we can use the FSS analysis. We should note that E p 
follows the FSS near the critical point, 

E p (p,L) ~ X{tL x l v ), t = (p c -p)/p c , (2.2) 

as far as the corrections to FSS arc negligible, where p c 
is the critical value of p for the infinite system (L — > oo) 
and v is the correlation-length critical exponent. Then, 
we can estimate p c from the size dependence of p c (L) 
using eq. (2.2) and, in turn, estimate T c through the 
relation p c = 1 — e - 2 J/ k BT c _ 

2.2 Percolating condition 

There are several choices of criterion to determine per- 
colating. For example, Machta et al. 14 ^ used both the 
extension rule and the topological rule for their stopping 
condition in the IC algorithm. The former rule is that 
some cluster has maximum extent L in at least one of the 
d directions in d-dimensional systems. The latter rule is 
that some cluster winds around the system in at least 
one of the d directions. We may use these percolating 
conditions in our PCC algorithm. Actually, we can use 
any rule to determine percolating, but FSS functions for 
physical quantities, therefore p c (L), depend on the rule. 

2.3 Distribution of p 

Let us consider the distribution of p, f(p). We change 
p based on the observation whether the system is per- 
colating or not, which leads to the negative feedback 
mechanism. Since the existence probability E p (p) is the 
probability that the systems percolates, the transition 
probabilities W from p to p + Ap and from p to p — Ap 



in the PCC algorithm are written as follows: 

W p ^ p+ a p = 1 - E p (p), 

W p ^ p -Ap = E p (p). 

In the vicinity of p c (L), we may employ the linear ap- 
proximation for E p (p), such as 

E p (p) = ^+a(p-p c (L)), (2.4) 

where a is the value of dE p /dp at p c {L). Then, this prob- 
lem is nothing but the Ehrenfest model for diffusion with 
a central force. 15 ' 16 ) In the steady state, the distribution 
of p, f(p), satisfies the relation 

f{p) = W p -A P ^ P f{p - Ap) + W p+A p^pf(p + Ap). 

(2.5) 

In the linear approximation, cq. (2.4), p takes the val- 
ues between p c {L) — l/2a and p c (L) + l/2a. Substi- 
tuting eqs. (2.3) and (2.4) into eq. (2.5) and denoting 
p = Pc(L) + iAp, we can show that f(p) is the binomial 
distribution, that is, 

f(p) cx n C n/2+i , (2.6) 

where n = (l/a)/Ap, i e [— n/2,n/2], and n C„/ 2 +i are 
the binomial coefficients. Thus, for large n, or small 
Ap, the distribution function f(p) becomes the Gaus- 
sian distribution with the average p c {L) and the variance 
a 2 = (n/4) (Ap) 2 — Ap/4a. For smaller Ap, the width 
of the distribution becomes narrower as a cx \/Ap. Since 
a is the value of dE p /dp at p c (L), we expect a cx L x l v 
using the FSS assumption, eq. (2.2). Thus, for larger L, 
the width of f(p) becomes narrower. 

In the Letter, 12 ) we have checked that the distribution 
of p actually approaches the Gaussian distribution with 
very narrow width for the 2D Ising model, and the re- 
sulting energy histogram is indistinguishable from that 
by the constant-temperature calculation. Moreover, we 
have obtained the expected Ap- and L-dependence for 
the width of f(p). 

2.4 Determination of next p 

In the process 3, we decrease or increase p by Ap. The 
difference Ap is a free parameter in our algorithm. In the 
limit of small Ap we approach the canonical ensemble, 
but it takes a long time to equilibrate for small Ap. Using 
the same approximation as in the previous subsection 
and assuming that the subsequent steps are independent, 
we can show that the deviation of the average value of p 
from p c (L) becomes smaller as a geometric progression 
with time. Since the geometric ratio is given by 1 — 2a Ap 
in this approximation, the convergence becomes slower 
for smaller Ap. Practically, we may start with rather 
large Ap, and gradually decrease Ap with monitoring 
the trail of the values of p. Small steps of preparation 
are enough for equilibration. 

We change p by every Monte Carlo step in our orig- 
inal proposal. 12 ^ As another way, we may measure the 
existence probability E p for a short time interval with 
keeping p constant, and then change p for the next short 
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time interval. In this process, recording the values of p 
and whether percolation occurred for each, we may ex- 
tract the information on the final p c {L) efficiently using 
the Bayesian statistics. 18 - 1 A more deterministic way of 
adjusting p c (L) may be considered. Solving E p (p) = 1/2 
iteratively by the Newton method may lead to the ad- 
justment of p c (L). 19 ) It is quite interesting to improve 
and extend the method of PCC algorithm. 

2. 5 Checking value of E p 

We have chosen the value of E p which gives p c {L) as 
1/2 because it is the simplest. In case the critical value 
of E p at the critical point of the infinite system is far 
from 1/2, it is convenient to employ the checking value 
of E p different from 1/2. We may modify the update 
process such that this value is different from 1/2. If we 
want to choose the checking value of E p as e p , we may 
modify the process 3 as follows: 
3'. If the system is percolating in the previous test, 
decrease p by Ap with the probability s = 
min[(l/2)/e p , 1] and increase p with the probabil- 
ity 1 — s. On the contrary, if the system is not 
percolating, increase p by Ap with the probability 
s' = min[(l/2)/(l — e p ), 1] and decrease p with the 
probability 1 — s'. 
In this way, we can control the checking value of E p . 

3. Results 

3.1 3D Ising model 

Here we present the result of the 3D Ising model. We 
have simulated the Ising model on the simple cubic lat- 
tice by using the PCC algorithm. We have treated the 
systems with linear sizes L=8, 12, 16, 24, 32, 48, and 64. 
For the criterion to determine percolating, we have em- 
ployed the topological rule 14 - 1 in the present study. After 
10,000 Monte Carlo sweeps of determining p c (L) with 
gradually reducing Ap, we have made 100,000 Monte 
Carlo sweeps to take the thermal average; we have made 
10 runs for each size to get better statistics and to 
evaluate the statistical errors. We have started with 
Ap = 10~ 3 for L = 8 and Ap = 10~ 4 for L = 64. For 
the intermediate sizes, we have started with Ap between 
these two values. The final value of Ap has been chosen 
as Ap = 10 -5 x L -1 - 6 for the system size L. Actually, 
the schedule of decreasing Ap is not so serious. 

We plot the size-dependent T C (L) as a function of l/L 
for the 3D Ising model in Fig. 1. From now on, we rep- 
resent the temperature in units of J/ks- The error bars 
are within the size of the mark. The critical tempera- 
ture T c can be estimated by the FSS relation, eq. (2.2). 
Including the corrections to FSS, we have 



T C (L) =T c + aL- 1 /»{l + bL-"), 



(3.1) 



where T C (L) is given through p c (L) = 1 - I^bT c {V) ^ 
and lo is the correction-to-FSS exponent. Since there 
are five fitting parameters in eq. (3.1), it is not easy to 
get accurate estimates of the critical point and critical 
exponents. In a high-resolution Monte Carlo study, Fer- 
renberg and Landau 17 ) employed a FSS analysis to get 




Fig. 1. Plot of T C (L) (in units of J/ks) as a function of l/L for 
the 3D Ising model. The system sizes are L = 8, 12, 16, 24, 32, 48, 
and 64. 



accurate estimates. They first determine the exponent 
v, and with v determined quite accurately they then es- 
timate T c . Since their procedure is well fitted for our 
algorithm, we use the same idea. 

We first note that the inverse-temperature derivatives 
of the logarithm for the moment of magnetization to obey 
the following FSS relations, 
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= a'L 1 ! v {l + b'L-"), 



(3.2) 



where K = 1/ksT and (•••) denotes the thermal av- 
erage. These FSS relations hold for the values of to at 
the size-dependent T C (L) as well as those at the fixed T c 
for the infinite system. Here, we treat the variables at 
T C {L). Since we calculate the left-hand side of eq. (3.2) 
by the general formula to calculate the if-derivative of 
any quantity {A), 



d(A) 
dK 



{A) (E) - (AE) 



(3.3) 



where E is the energy, we can extract v without de- 
termining T c for the infinite system. The inverse- 
temperature derivative of the Binder parameter, 20 ) 
which is defined as 



<" 4 ) 
(to 2 ) 2 



(3.4) 



also obeys the FSS relation as in eq. (3.2). What we do 
is that we measure the inverse-temperature derivatives 
of In (to™) and g at the size-dependent T C (L) and make 
an analysis based on the FSS relations, eq. (3.2). 

We plot the derivatives dg/dK, d\n(\m\) /dK and 
<91n (to 2 ) /dK at T C (L) as a function of L in logarithmic 
scale in Fig. 2. The error bars are again within the size of 
the mark. We find the power-law size dependence from 
the linearity of the data with small corrections to FSS; 
wc estimate the exponent 1/v = 1.594(8) using eq. (3.2). 
Here, the number in the parenthesis denotes the uncer- 
tainty in the last digit. Wc have used the average of three 
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Fig. 2. Plot of derivatives at T C (L) as a function of L for the 3D 
Ising model in logarithmic scale. 



Fig. 4. FSS plot of f m (m) for the 3D Ising model, where /3/u = 
0.517. The system sizes are L = 16, 32, and 64. The inset shows 
the raw data. 




Fig. 3. Plot of (\m\) at T C (L) as a function of L for the 3D Ising 
model in logarithmic scale. 



data. The correction-to-FSS exponent uu is 1.2(5), which 
is a little bit larger than the recent value u = 0.87(9). 21 ) 
Our estimate of T c using cq. (3.1) is T c = 4.5108(7), 
that is, 1/T C = 0.22169(4). Both estimates of T c and v 
are consistent with the estimates of the recent study, 17 ) 
l/v = 1.590(2) and 1/T C = 0.2216595(26). The solid 
curve in Fig. 1 is the best fitted curve for eq. (3.1) with 
v determined accurately first. 

In order to discuss another exponent, we plot the aver- 
age of the magnetization (|m|) at T C (L) as a function of 
L in logarithmic scale in Fig. 3. We use the FSS relation 
with the corrections to FSS, 



<H> 



T=T C (L) 



= cL-M"(l + dL- 



(3.5) 



for the estimate of the magnetization exponent (3. From 
the least square fit using eq. (3.5), we have (3/v = 
0.517(8), which is again consistent with the recent es- 
timate, 17 ) (3/v = 0.518(7). 

It is interesting to study the distribution function of 
physical quantities. We show the FSS plot of the dis- 



tribution function f m {m) in Fig. 4, based on the FSS 
relation, 



f m (m) T=Te(L) ~ L^h(mL^). 



(3.6) 



The inset of Fig. 4 shows the raw data of the distribution 
functions / m (m) for linear sizes L = 16,32, and 64. The 
scaled data show very good FSS behavior; that is, the 
data of different sizes are collapsed on a single curve. 

3. 2 3D bond percolation problem 

The idea of the PCC algorithm is based only on the 
property of a percolation problem. Thus, it is straight- 
forward, or even easier, to apply this algorithm to the 
geometric percolation problem. The partition function 
for the bond percolation problem is written as 



G'CG 



KG') 



(1-P) 



N b -b(G') 



(3.7) 



Here G is all the configurations, or the graph, and b(G') 
is the number of occupied bonds in the subgraph G' . 
The sum is over all subgraphs G' of G. The probability 
of bond occupation is denoted by p, and iVj, is the total 
number of bonds in the system. In the bond percolation 
problem, wc arc to locate the percolation threshold p c . 
We change p by the small amount of Ap in the process of 
simulation, and determine the size-dependent p c (L) au- 
tomatically as in the case of the Ising model. In the limit 
of Ap — > 0, it becomes the usual percolation problem. 
One thing we should have in mind is that we determine 
whether the bond is occupied or not with the probability 
p for each bond. 

We have studied the 3D bond-percolation model with 
linear sizes L = 8, 12, 16, 24, 32, 48, and 64. Almost the 
same conditions are used as in the 3D Ising model. We 
plot the size-dependent p c (L) clS cL function of 1/L in 
Fig. 5. The FSS relation for p c (L) is given by the equa- 
tion similar to eq. (3.1), but we follow the same scheme 
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Fig. 5. Plot of Pc(L) as a function of 1/L for the 3D bond per- 
colation problem. The system sizes are L = 8, 12, 16, 24, 32,48, 
and 64. 




as in the Ising model to get better estimates of the crit- 
ical point and exponents; that is, we first estimate the 
critical exponent v. The fraction of lattice sites in the 
largest cluster c plays a role of the order parameter. 
Thus, we consider the p-derivative of the moments of 
c. The derivative of any quantity (A) with respect to p 
is given by 



Fig. 6. Plot of derivatives at p c (L) as a function of L for the 3D 
bond percolation problem in logarithmic scale. 



d(A) (Ab)-(A)(b) 



dp 



p(l-p) 



(3i 



where b is the number of occupied bonds, and (• • • ) de- 
notes the sample average. Equation (3.8) corresponds to 
the general formula for the if -derivative, eq. (3.3). We 
may derive eq. (3.8) starting from the expression for the 
partition function, eq. (3.7), or using the correspondence 
based on the KF relation, p = 1 — pr 2J / i,iT . 

We have calculated the logarithmic derivatives of (c) 
and (c 2 ), where c is the fraction of lattice sites in the 
largest cluster. Since the moment ratio (c) 2 / (c 2 ) has the 
same FSS property as the Binder parameter, eq. (3.4), 
we also calculate the p-derivative of this moment ratio. 
We plot the derivatives <9((c) 2 / (c 2 ))/<9p, <91n (c) /dp and 
dln(c 2 ) I dp at p c (L) as a function of L in logarithmic 
scale in Fig. 6. We find the power-law size dependence 
from the linearity of the data with small corrections to 
FSS; we estimate the exponent 1/v = 1.12(5) from the 
slopes of lines. The correction-to-FSS exponent u is 
1.1(5), which is a little bit smaller than the recent es- 
timate for site percolation problem, uj = 1. 62(13). 2V> We 
can use the FSS form similar to eq. (3.1) for the esti- 
mate of p c ; our estimate using the value of v is p c = 
0.24881(3). Both estimates of p c and v are consistent 
with the estimates of the recent study, 22 ' 1/v = 1.12(3) 
and p c = 0.2488126(5). The solid curve in Fig. 5 is the 
best fitted curve for p c , which is similar to eq. (3.1), with 
v determined accurately first. 

We plot the fraction of lattice sites in the largest clus- 
ter (c) at p c {L) as a function of L in logarithmic scale 
in Fig. 7. Since (c) has the same FSS form as the mag- 
netization, we use eq. (3.5) for the estimate of (3/u. Us- 
ing the least square fit, we have (3/v — 0.474(5), which 
is again consistent with the estimate of recent study, 22 ' 




Fig. 7. Plot of (c) at p c (L) as a function of L for the 3D bond 
percolation problem in logarithmic scale. 



(3/v = 0.476(5). 

Finally, we show the FSS plot of the distribution func- 
tion / c (c) in Fig. 8. The raw data of / c (c) for linear sizes 
L = 16, 32, and 64 are given in the inset of Fig. 8. Again, 
in the percolation problem, the data show very good FSS 
behavior. 

4. Summary and Discussions 

To summarize, we have given a detailed description of 
the newly proposed PCC algorithm. We have applied the 
PCC algorithm to the study of the 3D Ising model and 
the 3D bond percolation problem. We have employed a 
refined analysis of FSS, which uses the same scheme as 
suggested by Ferrenberg and Landau. 17 ' Our results for 
the Ising model and the bond percolation problem are 
consistent with those of the previous works. 17 ' 22 ' It is 
to be noted that we make simulations at a single crit- 
ical point p c (L) for each system size. Thus, we need 
much less efforts compared with the usual procedure for 
making simulations at several different temperatures to 



6 



Yusukc TOMITA and Yutaka Okabe 



o L=64 
1 - 9 L=32 




' J 1 " " 2 

cL p/v 

Fig. 8. FSS plot of / c (c) for the 3D bond percolation problem, 
where (3/v = 0.474. The system sizes are L = 16, 32, and 64. 
The inset shows the raw data. 



extract the critical point and to estimate critical expo- 
nents. 

We have estimated v from the temperature derivative 
of the moments. There is an alternative way to esti- 
mate v without determining T c . Let us calculate the 
size-dependent T c in two ways. We may use different 
criteria for percolating condition; we may use different 
values of E p which gives p c {L). Then, the difference of 
two T c (L)'s follows the FSS relation as 

TW(L) - T^ 2) (L) = d'L- x l v (\ + b"L~ u ), (4.1) 

which also leads to the direct determination of v. 

There are several directions for the application of the 
PCC algorithm. We can use the PCC algorithm to any 
problem where the mapping to the cluster formalism ex- 
ists. It is straightforward to apply this algorithm to 
the diluted Ising (Potts) models. The PCC algorithm 
is quite useful for investigating the self-averaging prop- 
erties of random systems, where the distribution of T c 
due to randomness is essential. It is because we can 
determine the sample-dependent p c (L) quite easily. We 
have already applied the PCC algorithm to the 2D site- 
diluted Ising model, 23 > and have studied the crossover 
and self-averaging properties. It is also interesting to ex- 
tend the PCC algorithm to the problem of the vector 
order parameter. We have already succeeded in apply- 
ing the PCC algorithm to the classical XY model, 24 ) and 
have shown that the PCC algorithm is useful not only 
for the analysis of the second-order transition but also for 
that of the transition of the Kosterlitz-Thouless type. 

In the PCC algorithm, the cluster representation is 
used in two ways. First, we make a cluster flip as in 
the SW algorithm. 10 ) Second, we change p depending 
on the observation whether clusters are percolating or 
not. However, the percolating properties arc not essen- 



tial. We have used the FSS relation for E p , eq. (2.2), 
to determine the critical point. We may use quantities 
other than E p which have the similar FSS relation with 

a single scaling variable. We could generalize the PCC 
algorithm for a problem where the mapping to the clus- 
ter formalism does not exist. For example, we may study 
the systems with a frustration by the generalized scheme 
of the PCC algorithm. The application of the PCC al- 
gorithm to quantum spin systems is also an interesting 
subject, and now in progress. 
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